home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Cream of the Crop 26
/
Cream of the Crop 26.iso
/
os2
/
octa209s.zip
/
octave-2.09
/
libcruft
/
linpack
/
dgbfa.f
next >
Wrap
Text File
|
1996-07-19
|
5KB
|
175 lines
SUBROUTINE DGBFA(ABD,LDA,N,ML,MU,IPVT,INFO)
INTEGER LDA,N,ML,MU,IPVT(1),INFO
DOUBLE PRECISION ABD(LDA,1)
C
C DGBFA FACTORS A DOUBLE PRECISION BAND MATRIX BY ELIMINATION.
C
C DGBFA IS USUALLY CALLED BY DGBCO, BUT IT CAN BE CALLED
C DIRECTLY WITH A SAVING IN TIME IF RCOND IS NOT NEEDED.
C
C ON ENTRY
C
C ABD DOUBLE PRECISION(LDA, N)
C CONTAINS THE MATRIX IN BAND STORAGE. THE COLUMNS
C OF THE MATRIX ARE STORED IN THE COLUMNS OF ABD AND
C THE DIAGONALS OF THE MATRIX ARE STORED IN ROWS
C ML+1 THROUGH 2*ML+MU+1 OF ABD .
C SEE THE COMMENTS BELOW FOR DETAILS.
C
C LDA INTEGER
C THE LEADING DIMENSION OF THE ARRAY ABD .
C LDA MUST BE .GE. 2*ML + MU + 1 .
C
C N INTEGER
C THE ORDER OF THE ORIGINAL MATRIX.
C
C ML INTEGER
C NUMBER OF DIAGONALS BELOW THE MAIN DIAGONAL.
C 0 .LE. ML .LT. N .
C
C MU INTEGER
C NUMBER OF DIAGONALS ABOVE THE MAIN DIAGONAL.
C 0 .LE. MU .LT. N .
C MORE EFFICIENT IF ML .LE. MU .
C ON RETURN
C
C ABD AN UPPER TRIANGULAR MATRIX IN BAND STORAGE AND
C THE MULTIPLIERS WHICH WERE USED TO OBTAIN IT.
C THE FACTORIZATION CAN BE WRITTEN A = L*U WHERE
C L IS A PRODUCT OF PERMUTATION AND UNIT LOWER
C TRIANGULAR MATRICES AND U IS UPPER TRIANGULAR.
C
C IPVT INTEGER(N)
C AN INTEGER VECTOR OF PIVOT INDICES.
C
C INFO INTEGER
C = 0 NORMAL VALUE.
C = K IF U(K,K) .EQ. 0.0 . THIS IS NOT AN ERROR
C CONDITION FOR THIS SUBROUTINE, BUT IT DOES
C INDICATE THAT DGBSL WILL DIVIDE BY ZERO IF
C CALLED. USE RCOND IN DGBCO FOR A RELIABLE
C INDICATION OF SINGULARITY.
C
C BAND STORAGE
C
C IF A IS A BAND MATRIX, THE FOLLOWING PROGRAM SEGMENT
C WILL SET UP THE INPUT.
C
C ML = (BAND WIDTH BELOW THE DIAGONAL)
C MU = (BAND WIDTH ABOVE THE DIAGONAL)
C M = ML + MU + 1
C DO 20 J = 1, N
C I1 = MAX0(1, J-MU)
C I2 = MIN0(N, J+ML)
C DO 10 I = I1, I2
C K = I - J + M
C ABD(K,J) = A(I,J)
C 10 CONTINUE
C 20 CONTINUE
C
C THIS USES ROWS ML+1 THROUGH 2*ML+MU+1 OF ABD .
C IN ADDITION, THE FIRST ML ROWS IN ABD ARE USED FOR
C ELEMENTS GENERATED DURING THE TRIANGULARIZATION.
C THE TOTAL NUMBER OF ROWS NEEDED IN ABD IS 2*ML+MU+1 .
C THE ML+MU BY ML+MU UPPER LEFT TRIANGLE AND THE
C ML BY ML LOWER RIGHT TRIANGLE ARE NOT REFERENCED.
C
C LINPACK. THIS VERSION DATED 08/14/78 .
C CLEVE MOLER, UNIVERSITY OF NEW MEXICO, ARGONNE NATIONAL LAB.
C
C SUBROUTINES AND FUNCTIONS
C
C BLAS DAXPY,DSCAL,IDAMAX
C FORTRAN MAX0,MIN0
C
C INTERNAL VARIABLES
C
DOUBLE PRECISION T
INTEGER I,IDAMAX,I0,J,JU,JZ,J0,J1,K,KP1,L,LM,M,MM,NM1
C
C
M = ML + MU + 1
INFO = 0
C
C ZERO INITIAL FILL-IN COLUMNS
C
J0 = MU + 2
J1 = MIN0(N,M) - 1
IF (J1 .LT. J0) GO TO 30
DO 20 JZ = J0, J1
I0 = M + 1 - JZ
DO 10 I = I0, ML
ABD(I,JZ) = 0.0D0
10 CONTINUE
20 CONTINUE
30 CONTINUE
JZ = J1
JU = 0
C
C GAUSSIAN ELIMINATION WITH PARTIAL PIVOTING
C
NM1 = N - 1
IF (NM1 .LT. 1) GO TO 130
DO 120 K = 1, NM1
KP1 = K + 1
C
C ZERO NEXT FILL-IN COLUMN
C
JZ = JZ + 1
IF (JZ .GT. N) GO TO 50
IF (ML .LT. 1) GO TO 50
DO 40 I = 1, ML
ABD(I,JZ) = 0.0D0
40 CONTINUE
50 CONTINUE
C
C FIND L = PIVOT INDEX
C
LM = MIN0(ML,N-K)
L = IDAMAX(LM+1,ABD(M,K),1) + M - 1
IPVT(K) = L + K - M
C
C ZERO PIVOT IMPLIES THIS COLUMN ALREADY TRIANGULARIZED
C
IF (ABD(L,K) .EQ. 0.0D0) GO TO 100
C
C INTERCHANGE IF NECESSARY
C
IF (L .EQ. M) GO TO 60
T = ABD(L,K)
ABD(L,K) = ABD(M,K)
ABD(M,K) = T
60 CONTINUE
C
C COMPUTE MULTIPLIERS
C
T = -1.0D0/ABD(M,K)
CALL DSCAL(LM,T,ABD(M+1,K),1)
C
C ROW ELIMINATION WITH COLUMN INDEXING
C
JU = MIN0(MAX0(JU,MU+IPVT(K)),N)
MM = M
IF (JU .LT. KP1) GO TO 90
DO 80 J = KP1, JU
L = L - 1
MM = MM - 1
T = ABD(L,J)
IF (L .EQ. MM) GO TO 70
ABD(L,J) = ABD(MM,J)
ABD(MM,J) = T
70 CONTINUE
CALL DAXPY(LM,T,ABD(M+1,K),1,ABD(MM+1,J),1)
80 CONTINUE
90 CONTINUE
GO TO 110
100 CONTINUE
INFO = K
110 CONTINUE
120 CONTINUE
130 CONTINUE
IPVT(N) = N
IF (ABD(M,N) .EQ. 0.0D0) INFO = N
RETURN
END